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. Abstract. We study the evolution of an embedded protoplanet in a circumstellar disk using the 3D-Radiation Hydro 

(f) ■ code TRAMP, and treat the thermodynamics of the gas properly in three dimensions. The primary interest of 

this work lies in the demonstration and testing of the numerical method. We show how far numerical parameters 
can influence the simulations of gap opening. We study a standard reference model under various numerical 
approximations. Then we compare the commonly used locally isothermal approximation to the radiation hydro 
simulation using an equation for the internal energy. Models with different treatments of the mass accretion 
f*^ , process are compared. Often mass accumulates in the Roche lobe of the planet creating a hydrostatic atmosphere 

O i 1 around the planet. The gravitational torques induced by the spiral pattern of the disk onto the planet are not 

strongly affected in the average magnitude, but the short time scale fluctuations are stronger in the radiation 
hydro models. 

An interesting result of this work lies in the analysis of the temperature structure around the planet. The most 
striking effect of treating the thermodynamics properly is the formation of a hot pressure-supported bubble around 
the planet with a pressure scale height of H/R ~ 0.5 rather than a thin Keplerian circumplanetary accretion disk. 
' We also observe an outflow of gas above and below the planet during the gap opening phase. 
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1. Introduction 

During their formation process, massive, gaseous planets 
are believed to undergo a phase of evolution where they 
are still embedded in the protoplanetary disk. The gravi- 
tational influence of the planet onto the ambient disk leads 
to such features as spiral arms and, for planets sufficiently 
massive, an annular gap at the planetary radius. The back- 
reaction of the perturbed disk onto the planet generates 
torques, which induce a change of the orbital elements 
(semi-major axis and eccentricity) of the planet. The lin- 
ear phase for small planetary m asses has been studied 
for tw o-di mensional flat disks by iGoldreich fc Tremainel 
l)l98fj(l and IWardl ill 9971), and m ore recently in three di- 
mensions bv lTanaka et alJ l|2002|) . Analytical formulae for 
the (inward) radial migration rate have been obtained. 
These linear (local) analyses have been supplemented by 
a variety of numerical simulations that also consider the 



global disk and larger planetary masses. Tw o-dimensional 
flat d is ks wi t h plan ets w ere modeled by iBrvden et alJ 
l|l999j) : iKlevI <|l999j) and iLubow et all l|l999l) applying 
a locally isothermal assumption for the radial tempera- 
ture distribution. Later, global three-dimensional simula- 
tions of embedded pla nets were performe d on equidistant 
grids (Kiev et al.l200lt) . with nested grids I D'Angelo et al 



2003b ) and by applying local grid refinement ijBate et al 
20031) . They confirmed the linear results for very small 
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planets as well as the estimates for gap opening large 
planets. In these simulations the planet is held on a 
fixed circular orbit and the mass accretion and migration 
rate are estimated from the obtained density distribution. 
Simu lations with a moving single planet (|Nelson et alJ 
|2000|) are again in agreement with the inward migration 
scenario. 

These hydrodynamical simulations have been ex- 
tended recently to full 3D -MHD calculations though 
neglecting vertical gravity llPapaloizou fc Nelsonl l2003t 
iNelson fc Papaloizoul2003tlWinters et all2003[) . Here, the 
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turbulent state of the disk is modeled directly, and the mi- 
gration rate is found to fluctuate strongly around a mean 
value. 

All the above mentioned simulations treat the disk 
as isothermal loc ally, althou gh polytro pic disks have 
been considered in iKlevI l|l999|) . Recently. D'Aneelo et all 
<l2003al) added thermodynamical effects in a two- 
dimensional model. They assume that the dissipation pro- 
duced in the disk (by viscosity) is radiated locally while 
the equation of state includes dissociation and ionization 
effects. Using this method in addition to nested grid re- 
finement around the planet it also has been possible for 
the first time to obtain an estimate of the two-dimensional 
temperature structure in the disk's midplane in the vicin- 
ity of the planet. However, these simulations still suffer 
from their two-dimensionality. 

In the present paper we extend the radiation hydrody- 
namical planet-disk calculations to three dimensions, by 
applying a flux-limited diffusion approximation for the 
radiation. This allows us to study directly the dynamical 
influence of the planetary accretion luminosity and to de- 
termine the three-dimensional temperature structure in 
the vicinity the planet. 

Here, we focus on numerical issues, and study the ef- 
fects of the radiation transport, of resolution and mass ac- 
cretion. By introducing radiation and opacities the model 
is no longer scale-free, as it is in the isothermal case. 
Hence, there is a difference, whether the planet is located 
at 5 AU or 0.5 AU, say. In addition to influencing the 
gravitational torques acting on the planet, the disk mass 
now has an effect on the temperature structure as well, 
because it provides the optical depth. 

We find that for the test parameters (1 Jupiter mass 
at 5 AU) and the low numerical resolution that we apply, 
there is no significant acceleration or breaking of the mean 
radial drift in comparison to the isothermal approxima- 
tion. However, the fluctuations in the radial drift appear 
an order of magnitude larger in the radiative models with 
realistic accretion than in the isothermal models. 

We find the circumplanetary "accretion disk" to be 
rather thick, i.e. the pressure scale height is H/R k, 0.5 in 
our models, which is even thicker th an the results from 2D 
simulations <|D'Angelo et al The deviation from a 

Keplcrian profile in this "disk" is about 50%, which also 
indicates the strong support from the thermal pressure. 
Under such conditions the thin disk approximation breaks 
down and one has to deal with a circumplanetary cloud 
filling a sizable fraction of the Roche lobe rather than a 
disk. 

The detailed modeling of the observational appearance 
of a planet interacting with a disk and a detailed param- 
eter study are left to follow up studies. 

In the next section we describe our model setup, the 
boundary and initial conditions. In section 3 we present 
an analysis of our models, and we conclude in section 4. 



2. Model characterization 

We consider a Jupiter sized protoplanet embedded in a 
protoplanetary disk in three dimensions, where we use 
a spherical polar coordinate system (r,9,(p). The central 
star lies at the origin of the coordinate system and the 
disk midplane coincides with the 9 = plane. We work in 
a coordinate frame corotating with the planet, and use a 
special formulation avoiding explicit source terms to con- 
serve angular momentum locally and globally, see IKlevI 

(EI - 



2.1. Physical setup 

All models we consider in this paper have a radial ex- 
tent from 1.25AU to 20AU, where a Jupiter mass planet 
is located at 5 AU. The vertical extent is from —10° be- 
low the midplane to +10° above the midplane of the disk. 
In the azimuthal (^-direction we consider a complete an- 
nulus from to 2n. The initial temperature distribution 
is proportional to r _1 and "vertically" (i.e. in the polar 
direction) constant. This temperature is chosen to repro- 
duce a local pressure scale height H/r = 0.05 for all radii. 
The actual vertical extent of the computational domain 
is then approximately ±3.5 x H/r. The initial density 



distribution is radially proportional to 



with a 



plane density at 5 AU of p(5AU 7 9 = 0) = 10 -11 g cm" 3 . 
This corresponds approximately to a local surface density 
S = 2Hp = 75 [gem -2 ] at 5 AU, i.e. the l ocation of the 
planet . This i s simi l ar to the valu e used bv lBrvden et alJ 
(ll99^ : ll^lll99^ : lLubowet alJ lll999h . 

Due to the assumed local vertical isothermal disk as 
initial setup, the vertical density structure is given to first 
order by a Gaussian: 



p(r, 9) = p(r, 0) e 



(i) 



The vertical ve and radial v r velocities vanish initially. 
The initial azimuthal velocity is given by the equilibrium 
of gravity, centrifugal acceleration and the radial pres- 
sure gradient. This is not an exact solution because the 
Gaussian density distribution is only valid in the thin disk 
approximation. Thus, as we resolve the disk out to ±10° 
a small discrepancy appears. Nevertheless, the model is 
only slightly out of equilibrium and will relax quickly into 
a more stable state, which deviates only marginally from 
the initial state. 

In the models we do not use any explicit physical disk 
viscosity and model the ideal hydrodynamical equations. 
However, we use an artificia l bulk visco sity in the form of 
lyon Neumann & Richtmverl l|l95fl() . see IStone fc Normanl 
(|l992l) . Applying an explicit viscosity, e.g. a-viscosity 
would have introduced one further free parameter, which 
we seek to avoid. Although MHD turb ulence is a good can- 
didat e to produce an a-like viscosity ijShakura fc Sunvaevl 
Il973h . it is not obvious that a would be a constant in time 
and space (see e.g. IWimsch et al.ll2005|) . It is also unclear 
in how far the magneto-rotational instability would oper- 
ate in the gap region and in the circumplanetary bubble. 
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A second reason for not using an explicit viscosity is 
that it would have led to additional heating on top of the 
adiabatic contraction of the circumplanetary bubble, thus 
making it more difficult to identify the importance of the 
heating by the accretion onto the planet. 

2.2. Gravitational potential 

The Roche lobe around the planet determines the region 
where mass is gravitationally bound to the planet. The 
Roche lobe is pear-shaped with a sharp edge at the side 
towards the star and thus is usually approximated by the 
Hill sphere. This sphere is centered at the location of the 
planet and completely includes the Roche lobe. We will 
usually refer to the Roche lobe and its approximate size 
given by the Hill radius. 

To calculate the potential due to a planet of mass M p 
at the radius r p in the midplane of the accretion disk, 
we determine the radius ah of the Hill sphere around 
the ^planet, which is a function of the reduced mass fi 

M Q +M % 



of the system 



a h 



fl\ 3 



(2) 



For distances from the planet larger than a critical fraction 
q g of the Hill radius ah we use the exact gravitational 
potential <I> as the sum of stellar $ Q and planetary $ p 
potential 



<f> = 5v 



% = 



M P G 



(3) 



where r is the radius vector pointing from the star to a 
specific location in the disk. For closer distances we use 
a cubic expansion which is 2 nd order accurate at a critical 
distance a g = g g ah (i.e. it agrees with the outer potential 
in its 1 st and 2 nd derivative) and then becomes shallower 
towards the center with a vanishing gradient at the planet 
location: 



-M p G 



- r v) 3 o - r P y 



2 



(4) 



The benefit of such a treatment for the smoothing is that 
the potential is accurate down to the critical distance, 
whereas using the more classical approach 



M P G 



\J(r- r p f + a 



(5) 



leads to a significant underestimation of the potential 
depth already at a g . 

In our models we test the influence the parameter q g 
has on the mass of gas that accumulates in the Roche 
lobe of the planet. 

2.3. Accretion 

Accretion is treated as in iKlevI <ll999tK where a certain 
fraction of mass is taken out per time step within a given 



accretion radius a acc = a acc ah around the planet. Here, we 
take this fraction q acc of the Hill radius as a free param- 
eter and study its influence on the gas density and tem- 
perature in the vicinity of the planet. We believe that our 
models that treat accretion and conserve the total energy 
(DR and DR4) have the highest relevance for the tempera- 
ture structure of the circumplanetary cloud. These models 
yield a lower estimate for the temperatures, because using 
a deeper potential might still increase the temperatures. 
We also measure the mass contained in the Roche lobe and 
the mass change in the Roche lobe per time unit, which 
is a good tool to show that the planet is accreting at a 
steady rate and that the system is in equilibrium, i.e. no 
mass is piling up in the Roche lobe any longer. 

2.4. Numerical scheme 

We use the TR AMP hydro code as described in 
iKlahr et all l)l999|) . TRAMP is a versatile 3D radiation 
hydro package that can use various advection schemes. 
For the purp ose of this paper we apply a monotonic trans- 
port scheme ((van Leerlll977j) . The radiation transport in 
TRAMP is treated in the flux- limited diffusion approxi- 
matio n using th e flux li miter s of lLevermore fe Pomranind 
(Il98ll) . see also IKlevI df989h. and dust opacities typical 
for the solar nebula ((Bell fc Linl The scheme to 

tr eat internal e nergy and radiation energy density is as 
in iKlev fc Linl l(l99rl) . Adding the equations for internal 
and radiation energy and assuming that the latter is much 
smaller than the internal energy, which is certainly a very 
good approximation for protostellar disks we arrive at the 
following evolution equation for the internal energy 



d t e + V • (ue) = -pV • u + $ - V • F 



(0) 



where u is the velocity vector, p the thermal pressure, 
$ is the dissipation by viscosity. Here we use only artifi- 
cial viscosity in shocks to achieve the proper entropy gen- 
eration but no explicit physical viscosity in the smooth 
flow. Nevertheless, tensor viscosity is fully implemented 
in TRAMP and we can use it in future simulations of 
gap opening. For this paper containing test simulations 
we avoid the additional effects of viscous heating in or- 
der to study more carefully the accrctional heating in the 
vicinity of the planet. 

The radiation flux F is calculated in flux-limited dif- 
fusion approximation 



F = ^V£ R 
pn 



(7) 



with the Rosseland mean opacity k and the speed 
of light c. The flux lim i ter A is defined according to 
iLevermore fc Pomranind l|l98l|) . Here we assume local 
thermal equilibrium and set 



E B . = aT 4 . 



(8) 
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Numerically, the radiative diffusion part is solved sepa- 
rately as part of an operator splitting technique. The dif- 
fusion equation for 

d t E R = -VF = V — V£r (9) 

pK 

is solved implicitly to avoid time step limitations using e.g. 
the standard SOR (successive over-relaxation) scheme. 
The new temperature is then calculated from Eq. (jSJ. 

The code ha s been extensively tested in the context of 
accretion disks ((Klahr et al.lll993) . The new part that is 
implemented to the code is the inclusion and modification 
of the gravitational potential of the planet as described 
above. 

2.5. Boundary conditions 

The boundary conditions are chosen to minimize the ef- 
fects of wave reflections at the boundaries. In the polar 
and radial direction these are outflow conditions not al- 
lowing for inflow. Thus, for pressure, density and orthog- 
onal velocities, we apply a zero gradient condition, while 
the normal velocity components have a vanishing gradient 
only if the direction of the flow is outward of the com- 
putational domain. For inward directed flow the normal 
velocities components are set to zero. For the angular ve- 
locity we assume a Keplerian flow. No sponge layer as in 
iKlahr fc Bodenheimerl (^2003^ is used. In the azimuthal di- 
rection we apply periodic boundary conditions. 

3. The models 

The purpose of all the different models using the same 
physical parameters, i.e. stellar, planetary and disk mass 
as well as planet location, is to study the influence of the 
numerical treatment of the problem and to find the nu- 
merical parameters suited to study the gap opening as 
realistically as possible. 

In Table Q we give an overview of all the models we 
compare in the following. The second column (grid) states 
the numerical resolution (n r ,n^,n v ) of the grid. There 
are low resolution (100x20x200) models with a logarith- 
mic radial spacing and equidistant in the azimuthal and 
vertical direction. We also calculated some high resolution 
models, where the grid is in all dimensions logarithmically 
concentrated (indicated by Ic) towards the location of the 
planet, which results in a resolution in the Roche lobe up 
to four times higher than in the low resolution case. In 
these models there are about 2400 grid points within the 
Roche lobe. We tested the influence of the logarithmically 
distorted high resolution grid on a disk without a perturb- 
ing planet. Despite the grid refinement around the planet 
and quite distorted grid cells towards the edges of the grid, 
there are no unstable growing modes. 

Some models are calculated in the standard locally 
isothermal approximation (EOS = iso), some treat the 
internal energy equation plus flux limited radiation trans- 
port (EOS = rad). We also state the parameters for the 



smoothing of the gravitational potential q g and for the 
mass accretion onto the planet q aC c- Some models start 
with no gap initially, others use an axisymmctric initial 
model with a gap (ini = G). 

After briefly mentioning a test model for the non- 
equidistant high resolution grid, we will go through our 
models in the following order: first isothermal models, then 
radiative models, followed by non-accreting models, and 
finally accreting and energy conserving models. 




10 10 
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Fig. 1. Model CI: Snapshots taken every 26 orbits from 
t=60 to t=320 orbits, which is after the convergence of 
model CI. The radial slices of the density, temperature 
and azimuthal velocity are taken in the midplane at the 
location of the planet. Model CI defines the initial values 
for all the models that start off with an axisymmctric gap. 



3.1. Isothermal Models 

3.1.1. The axisymmetric model CI 

This model uses only one grid-cell in the azimuthal direc- 
tion but it includes the full potential of the planet, and 
mass is taken out of the grid from the center of the Roche 
lobe. It is used to run a model quickly into an initial state 
with a gap, which forms here not due to the gravitational 
torques by the planet but rather by the perturbation in 
gravity as it would result from a massive ring. However, 
as we shall see the obtained gap size and shape are in 
reasonable agreement with the realistic 3D case. 

In Fig.^we introduce the kind of plot that we will sub- 
sequently use for some of the following models. It shows 
quantities measured in the midplane at the azimuthal po- 
sition of the planet. Thus, it is one-dimensional data. We 
plot 10 snapshots, which are taken at equal time inter- 
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Fig. 2. Model CDI: Snapshots taken every 2 orbits from 
t=164 to t=184 orbits. The radial slices of the density are 
taken in the midplanc at the location of the planet. The 
surface density is azimuthally averaged. 
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Fig. 3. Model CDI: (upper row) accreted mass m a and 
accretion rate onto the planet rh a over time; (lower row) 
evolution of the mass contained in the Roche lobe ttiro 
and mass change in the Roche lobe mR, in units of Jupiter 
masses and orbital periods of the planet. 



vals, e.g. in the case of Fig.^every 26 orbits from t=60 to 
t=320 orbits. First we plot the density p. One clearly sees 
how the gap forms and how stable it is. Most of the disk 
is still at the initial values, which are given by the dashed 
line. The next plot shows the azimuthally averaged surface 
density S. This value is slightly higher than 75[gcm~ 2 ] at 
5 AU, because here we integrate over the total height of 
the computational domain and not only over ±iJ. The 
temperature T in the third plot does not fluctuate as this 
run is isothermal. 

The deviation of the local azimuthal velocity (v^/vk) 
from the Keplerian value is significant only in the vicin- 
ity of the "planet" . The shape of the rotational profile 
is determined by the gradient in the surface density and 
thus has no physical relevance for the real gap case where 
the gravitational torques are determining the shape of the 
gap. We also have to mention that this artificial gap struc- 
ture is not unique, but there is a family of solutions, where 
either the surface density profile or the azimuthal veloc- 
ity profile is a free parameter. Still this method generates 
a gap structure that is not too unrealistic to serve as an 
initial model for our simulations. 

The stationary state of model CI is given by the bal- 
ance of centrifugal forces, radial pressure gradients and 
the attraction of the combined stellar and ring potential, 
is used as the initial setup in all of the subsequent mod- 
els that include a C in their names (see Tab. Models 
without a C are started from the same initial condition as 
model CI is started, which is a disk without a gap. 




time [Orbits] 



Fig. 4. Model CDI: Evolution of planet migration rate in 
terms of migration time. The solid curve is based on the 
torques excluding the Roche lobe and the solid horizontal 
line is the time average. In addition we plot the torques 
based on an exclusion of 1.5 and 0.5 Hill radii (dotted and 
dashed lines respectively). In this particular model the gap 
is very clean, such that all three torques mostly coincide. 

3.1.2. Model CDI 

This is the standard isothermal three-dimensional model. 
Mass is taken out of the Roche lobe . Thus, the mode l 
is directly comparable to the ones in iKlev et aD l|200lh . 
The radial midplane profiles (see Fig. |5J) are measured at 
equal times between 164 and 184 orbits. The profile only 
changes very slowly. In the inner part of the disk one recog- 
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nizcs a decrease in density, which is related to two effects: 
First, the accretion process triggered by the tidal forces of 
the planet, which 'drives' the material through the inner 
boundary simulating accretion onto the star. This is pos- 
sible because we use non-reflecting boundary conditions 
at the radial boundaries, which allow for outflow but not 
for inflow. As an effect also the most outer grid cells loose 
a lot of mass, which shows at most in the high resolution 
models, e.g. model CD4I. 

The second reason for the mass depletion is the accre- 
tion process onto the planet, which occurs from both sides, 
from the inner as well as from the outer part of the disk. 
But the mass reservoir in the outer part of the disk is much 
larger, thus the decrease in mass is not visible over the 
sp an of the simu l ation. I n com p arison to th e simulations 
bv iBrvden et all l|l999j) . iKlevI l)l999j) . and iLubow et ail 
l)l999h the gap in the surface density is shallower. Those 
simulations used explicit tensor viscosity and started al- 
ready with a wider gap than we used. So in our case it will 
take many hundreds of orbits before the gap will reach a 
comparable depth. In terms of computational time such 
a run would take several months of computation time on 
our fastest available workstations. The code in its present 
version is not parallelized, due to the implicit solver of the 
radiation part. In the future we will consider runs with a 
much longer evolution time. 

Plotting the mass accreted onto the planet and the 
mass accumulated in the Roche lobe (see Fig.EJ) also shows 
that we are only slowly approaching a steady state of the 
solution. There is a mass of 2.2 x 10~ 5 Jupiter masses in 
the Roche lobe and the accretion rate, i.e. the mass taken 
off the grid, is 1.5 x 1CP 4 Jupiter masses per orbit. The 
time that mass resides in the Roche lobe in this model is 
only about a tenth of an orbital period. The non-zero ac- 
cretion rate is interesting because we include no viscosity 
in our model. One might expect that planets opening a 
gap in an inviscid disk could not accrete any more mass. 
There is the possibility that the planet-disk torques to- 
gether with the induced Reynolds stresses provide enough 
angular momentum transport to feed the planet via the 
circumplanetary accretion disk. Nevertheless, to clarify 
this issue of the accretion rate onto the planet the sim- 
ulations will have to be continued for many more orbits 
at high resolution. It may still be possible that accre- 
tion will stop after a thousand orbits, which is the case in 
non-viscous 2D simulations of gap opening. 

We also calculate the gravitational torques acting on 
the planet an d translate them into the typical radial mi- 
gration time l(Klev et al.ll200lh . We measure the torques 
in three ways: excluding the mass contained in the inner 
half Hill radius of the Roche lobe (in the following the 
dashed curves), excluding all the Roche lobe (solid curves) 
and excluding 1 5 times the Hill radius (dotted curves). 
iBate et all l|2003l) have used this method to give error bars 
to the measurements of the torques from the numerical 
simulations. This procedure is also very useful to study the 
relative importance of the torques generated from gas close 
to the planet. In Fig. 0] we see that the torques vary over 
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Fig. 5. Model CD4I: Snapshots taken each orbit from 
t=71 to t=81 orbits. The radial slices of the density are 
taken in the midplanc at the location of the planet. The 
surface density is azimuthally averaged. The dotted line 
indicates the values for model CDI for comparison. 

time but the mean drift rate corresponds to the 10 5 yrs ob- 
served in other simulations and which fit th e results from 
linear theory JGoldreich fc Tremaind Il980t IWardl Il997t 
iTanaka et alJl2002|) . It is remarkable how close the three 
measurements of the torques are at each time. This is a 
consequence of the fact that we take the mass out of a 
significant part of the Roche lobe producing a very clean 
gap. There is not enough matter left, to exert a torque on 
the planet. This changes in the models that do not treat 
mass accretion onto the planet (g a cc = 0.0) or remove gas 
from a smaller part of the Roche lobe (q acc = 0.2). 

3.1.3. Model CD4I 

We also have performed simulations with the same pa- 
rameters as for model CDI but with a higher resolution 
due to logarithmic concentration of grid cells around the 
planet. Now the gravitational smoothing length is set to 
a value of only 20% of the Hill radius (q g = 0.2), and the 
accretion radius to the same small region (g acc = 0.2). The 
radial slices (all orbits from t=71 to t=81 orbits) through 
the planet (see Fig. |SJ show that more mass has accumu- 
lated in the Roche lobe than in model CDI (the values 
from model CDI are given by the dotted line) . This result 
can be expected because of the deeper gravitational po- 
tential of the planet; the increase of density towards the 
planet is clearly visible. The inner disk is less depleted 
than in model CDI, but this indicates the earlier evolu- 
tionary step (81 orbits for model CD4I vs. 184 orbits for 
model CDI) . This model apparently loses more mass at the 
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outer boundary, but this is mainly due to the very coarse 
grid in combination with the open boundary conditions 
there. 

Nevertheless, the overall structure of the gap has not 
changed substantially. This model takes about 5 times 
more computational time than the low resolution case. 
This is mainly due to the fine resolution in the azimuthal 
direction, i.e. the Courant condition for the azimuthal ve- 
locity component. Table 1 shows that there is a little more 
mass accumulating in the Roche lobe (which was empty 
initially) than in the low resolution case. The torques are 
now strongly fluctuating while the average values are very 
similar to those of model CDI (see Table 1). The fluctua- 
tions are so strong that they even change their sign. Again 
the time averaged torques are insensitive to the radius 
around the planet, which has been excluded from the cal- 
culation of the torque. We suspect that the oscillations 
might be an effect of the high density concentration near 
the planet (see Fig.[5J), because small density asymmetries 
induce large torque variation. Similar fin dings have been 
obtain ed for nested grid simulations by iD'Aneelo et alJ 
ll20f)2h . 

In future work we will return to mesh refinement 
around the planet to study these effects in more detail. 

3.2. Radiative models 
3.2.1. Model CD 

The isothermal approach used in the previous models is 
a strong simplification of the thermodynamics of the disk 
near the planet and may not be well suited for studying 
the gap opening problem. For the radiative model CD we 
keep all parameters identical to model CDI, but use the 
equation for internal energy including flux-limited radi- 
ation transport, see Eq. |SJ Again, mass is removed from 
within the Roche lobe, which means that energy is taken 
out without releasing its accretion energy (this will be 
treated differently in models DR and DR4). The temper- 
ature and density structure of this model are displayed 
together with other models in Fig. [S] with dotted lines. In 
contrast to model CDI the disk has cooled dramatically 
because there is no internal heating in our inviscid model, 
while the region around the planet has become slightly 
hotter than in the isothermal gap case (Fig. EH). The rest 
of the gap is much cooler. There is also three times more 
mass piling up in the Roche lobe than in Model CDI (see 
also Table 1). This is a direct consequence of the radiative 
cooling. If the matter within the Roche lobe cools down 
it will contract leading to a higher density. The mass ac- 
cretion rate onto the planet is roughly three to four times 
larger than for the corresponding isothermal model. Thus, 
there is more mass stored in the Roche lobe, which also in- 
fluences the torque (see Fig- El - While the overall average 
migration time is only 10% less than in model CDI, there 
are now variations between the torques measured by ex- 
cluding more or less mass from the Roche lobe. Excluding 
cither the entire Roche lobe or even 1.5 times the Hill ra- 
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Fig. 6. Model CD4: Snapshot at t=5 orbits after the re- 
finement of the grid. Radial slices of density and temper- 
ature are taken in the midplane at the location of the 
planet, while surface density is the azimuthal average, 
(solid line = Model CD4; dotted line = Model CD; dashed 
line = Model CDS). 

dius gives very similar results, while excluding only 0.5 of 
the Hill radius doubles the amplitude of the fluctuations. 

This is the first hint that radiative cooling changes the 
structure of matter in the Roche lobe and consequently 
the resulting torques on the planet. 

3.2.2. Model CDS 

The following model is identical to the previous one but 
this time mass is taken out of only 20% of the Hill radius. 
Thus we can study the effect of the arbitrary size of this 
mass removal zone. As a result more mass accumulates 
in the Roche lobe close to the planet and the tempera- 
ture around the planet is even higher (Fig. dashed line). 
Nevertheless, the model is already strongly limited by res- 
olution as 20 % of the Hill radius means only two grid 
cells. 

The accretion rate is also close to the one in model CD. 
The average torques are still the same as in the isothermal 
case, model CDI. It also shows that the migration times 
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Fig. 7. Model CD: Evolution of the planet migration rate 
in terms of migration time. The solid curve is based on the 
torques excluding the Roche lobe and the solid horizontal 
line is the time average. In addition we plot the torques 
based on an exclusion of 0.5 Hill radii (dotted lines) and 
1.5 Hill radii (dashed lines). The horizontal lines refer to 
the time averaged values. 



are the longest if one excludes only half the Hill radius. 
This indicates the importance of the proper treatment of 
mass in the Roche lobe and its potential capability of al- 
tering the radial drift rate. This model is also used for 
comparison to model CD4 and CD in Fig. 

3.2.3. Model CD4 

Here, we increase the resolution and decrease the grav- 
itational smoothing again to 20% of the Hill radius (cf. 
model CD4I) . The radius out of which mass may be taken 
is decreased in the same fashion. Otherwise this model 
(see Fig. EJ is identical to models CD (dotted line) and 
CDS (dashed line). 

The most striking effect of the higher resolution is the 
higher temperature at the planet's location. Temperatures 
three times higher than in model CD are reached. 

This can be explained by two effects which are: 1. A 
deeper potential well leading to stronger compression and 
release of accretion energy. 2. The much higher densities 
near the planet leading to higher optical depth and lower 
cooling rates. 

Again the accretion rate is similar to model CD and 
the average torques are also not very different (see Table 
1). 

3.3. Non-accreting models 

In all of the previous models mass has been taken out 
from the grid in the vicinity of the planet. In the isother- 
mal models this might be a tolerable simplification as the 
average torques are not different between the high and low 
resolution cases. 



However, in the models with radiation transport one 
introduces a large error if one simply takes mass off the 
grid: Internal energy is removed from the simulation along 
with the mass, but the accretion energy released during 
final accretion phase onto the planet is neglected. One can 
now cither construct a model for this process and re-adjust 
the energy content of the gas after some mass has been 
taken away (see Models DR and DR4) or simply refrain 
from accreting mass and study the effects of the mass 
piling up. The first approach has one more free parameter 
than the second, which is the relative accretion rate (i.e. 
the fraction of the mass removed from the grid per orbit). 

We will first study more closely the effect of taking 
mass arbitrarily from the grid. Therefore, we go back to 
the isothermal model CDI but do not remove any mass 
this time. 

3.3.1. Model CDNI 

In comparison to model CDI the lack of accretion leads to 
a slightly higher concentration of mass around the planet, 
similar to the models with a deeper potential (see Table 
1). The mass in the Roche lobe is now three times greater 
than in model CDI. 

The major effect of not taking mass off the grid and 
letting it pile up in the Roche lobe is the introduction of 
strong fluctuations in the migration rate. The torques cal- 
culated including half of the Hill radius especially showed 
large fluctuations. Nevertheless, the time averaged migra- 
tion rate is still dominated by the structures outside the 
Roche lobe and thus are of the same magnitude as in all 
the previous models (see Table 1). 

3.3.2. Model CDN and CDN4 

In these models we start with a gap, do not remove any 
mass and perform the full radiation step. Models starting 
with a gap eventually converge to models without an ini- 
tial gap. This convergence can be demonstrated for the low 
resolution models (see Table 1), while the high resolution 
models are too expensive in computational effort to show 
this in a reasonable computation time. The general trend 
can be seen in Table 1. Density and temperature increase 
with the resolution of the Roche lobe while the torques are 
not much different in the high and low resolution cases. 

3.3.3. Model DN and DN4 

Instead of beginning with an artificial gap and then slowly 
filling the Roche lobe to accumulate mass around the 
planet, we start here with an unperturbed axisymmetric 
disk without a gap. The main result is that now a pres- 
sure supported atmosphere builds up around the planet, 
which becomes very hot. Both models become stationary 
and we observe maximum temperatures of about 600 K for 
model DN and 1150 K for the high resolution case DN4 
(see Table 1). The radial profiles for these models along 
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Fig. 8. Model DR4: Early stage of gap opening. Temperature in the r-9 plane of the protoplanetary disk at the 
azimuthal location of the planet. This snapshot is taken after 12 orbits, thus the gap is not completely empty yet. In 
this early phase the heating by the accretion onto the planet drives a strong convectional flow in the still optically 
thick disk material above the planet. Brightness is logarithmic temperature (lighter = warmer) ranging from 38Kto 
1060K and vectors give the logarithmic mass flux. We also plot iso-density lines with labels giving values in g cm~ 3 . 
Density ranges from 10 15 g cm~ 3 up to 1.3 x 10 10 g cm -3 . See the online edition for a color version of this plot. 



with the models DR and DR4 arc displayed in Fig. |^1 We 
remove no mass from the grid but treat the full radia- 
tive dynamics. The fact that the temperature is constant 
can be interpreted with a self regulated accretion process 
into the Roche lobe. The mass accretion rate balances the 
contraction of the cooling planetary envelope and this con- 
traction releases potential energy in the planetary gravi- 
tational field. Thus we have a balance between mass flux 
and radiation flux over the surface of the Roche lobe. The 
equilibrium between cooling and accretion now yields a 
temperature that only depends on the depth of the poten- 
tial well. This leads us to assume that the temperatures 
will be even higher if we can avoid the artificial smoothing 
of the potential well. This is partly tested by model DN4 
with higher resolution and a smaller smoothing parameter 
(see Table 1). Interestingly there is also a good agreement 
in temperature structure with the models that conserve 
total energy during accretion (see Fig.EJ. 

3.4. Accreting and energy conserving models 
3.4.1. Model DR and DR4 

This model is similar to model CD in the sense that we ac- 
crete mass from the Roche lobe. But this time it is started 
without a gap and we conserve the total energy locally 



in the grid cells where the gas is removed from the grid. 
Hence, the kinetic energy and the potential energy that 
can be released when the gas is accreted onto the planet 
are added locally to the internal energy. The radius of such 
a young planet probably will be about two times the cur- 
rent size of Jupiter (Bodenheimer, priv. communication). 

First we show the flow situation in model DR4 after 12 
orbits past the initial set up. The gap has not cleared com- 
pletely and the heating by accretion onto the planet drives 
strong vertical convection in the optically thick disk ma- 
terial above the planet (see Fig. EJ. This "fountain" flow 
is subsonic and will not escape from the disk because the 
velocity is too low in comparison to the escape velocity. 
Nevertheless, the effect might be very important for obser- 
vations because even embedded planets still too small to 
open a gap will probably produce this effect, which shows 
as an extended bubble of hot gas to the planet thus raising 
the local pressure scale height in the disk. 

Model DR approaches a quasi steady state after 100 
orbits (see Fig. |5J), which can also be inferred from the 
almost constant Roche lobe mass of 5 x 10 _4 A / /j up and the 
accretion rate onto the planet of 6 x lQ~ 4 Mj up per orbit 
(see Fie. ll0|) . Nevertheless, there are also short time scale 
fluctuations in the Roche lobe mass, which again show 
up in the fluctuations of the torques (see Fig. ITT|) . Based 
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Fig. 10. Model DR: (upper row) accreted mass m a and 
accretion rate m a onto the planet over time; (lower row) 
evolution of the mass contained in the Roche lobe m^ 
and mass change in the Roche lobe ttt-ro in units of Jupiter 
masses and orbital periods of the planet. 



10 



"[AU] 



Fig. 9. Model DR: Snapshot taken after t=121 orbits^ 
(solid line). Model DR4: after t=141 orbits (dotted line), u 
Model DN: after t=121 orbits (dashed line). Model DN4:^ 
Snapshot taken after t=55 orbits (dash dotted line).g 
Radial slices of density and temperature are taken in the $ 
midplane at the location of the planet, while surface den- J 
sity is the azimuthal average. ^ 
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on non-viscous 2D runs we know that it will take up to a 
thousand orbits before the gap is so deep that no accretion 
onto the planet occurs any longer. In 3D with radiation 
transport wc presently cannot afford such long runs. The 
resulting maximal temperature in the Roche lobe is 800 K 
for model DR and 1500 K for model DR4. A finer grid 
resolution might even lead to higher temperatures. 

For model DR we also display the density structure 
that results in the midplane (see Fig. I12|) after 122 orbits. 
The gap has clearly formed but still mass is left in the 
horseshoe orbit. This mass will decrease gradually over 
time. The density structure does not change significantly 
from any isothermal simulation, but the temperature of 
course does (see Fig. . The planet is the hottest spot 
in the simulation. We adopted a grey approximation for 
the radiation transport, thus continuum (wavelength de- 
pendent) radiation transport simulations (using our tem- 
perature and density structure) will show whether this hot 
circumplanetary gas-dust mixture will be observable. The 
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Fig. 11. Model DR: Evolution of planet migration rate in 
terms of migration time. The solid curve is based on the 
torques excluding the Roche lobe and the solid horizontal 
line is the time average. In addition we plot the torques 
based on an exclusion of 0.5 Hill radii (dotted lines) and 
1.5 Hill radii (dashed lines), plus the time averaged values. 



extent of the warm gas around the planet as well as the 
two mass feeding streams of warm gas hitting the circum- 
planetary cloud are clearly visible now. We do not resolve 
any circumplanetary accretion disk yet, for two reasons: 
First our resolution is too low and second the pressure 
dominates the hot gas around the planet strongly enough 
so that no Keplerian disk forms, and one is left with a 
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Fig. 12. Model DR: Density in the midplane of the 
protoplanetary disk after 121 orbits. Dark means low, 
light means high density. The density ranges from 2.4 x 
lCT 13 gcm~ 3 to 8.3 x 10 -11 gcm~ 3 . See the online edition 
for a color version of this plot. 



pressure-supported planetary envelope filling most of the 
Roche lobe. 

In Fig. El we display the pressure scale height of 
the gas in the midplane around the planet. In the cir- 
cumplanetary "disk" H p /R p is the inverse Mach number 
1/M — c s /vk of the flow around the planet. Due to the 
high temperatures the relative vertical thickness is about 
0.5, which is larger than those values typically used in thin 
Keplerian accretion disks with H/R < 0.1. In the case of 
H p /Rp = 0.5 this structure is better described by a rotat- 
ing envelope which additionally shows strong deviations 
from a Keplerian rotation law of about 50%. 

For model DR4 we display some details of the flow 
onto the planet (see Fig. I15f) . We plot the temperature in 
the r — 8 plane of the protoplanetary disk at the azimuthal 
location of the planet. The temperature is hottest around 
the planet and the gas in the gap is also heated by the 
accretion luminosity of the planet. The flow indicates that 
accretion mainly occurs via the poles of the planet and 
that there is convection in the envelope around the planet. 
Wc can see no inflow along the equatorial plane. The mass 
falling onto the central region around the planet originates 
only in part from the gap material above the planet. The 
largest contribution in mass is via the convective flow in 
the upper part of the circumplanetary cloud. 

This simulation also suggests that at least at this stage 
of gap opening and planet disk interaction there is no for- 
mation of moons from the circumplanetary material pos- 
sible. The reason is simple: with H p /R p = 0.5 and the 
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Fig. 13. Model DR: Temperature in the midplane of 
the protoplanetary disk after 121 orbits. Brightness corre- 
sponds to the logarithm of temperature. Black indicates 
10K and white 800K. See the online edition for a color 
version of this plot. 

resulting strongly sub-Keplerian rotation, all solid mate- 
rial will quickly rain down onto the planet. The maxi- 
mum radial drift velocity is given by the difference be- 
tween the actual azimuthal v elocity and and the Keplerian 
value l| Weidenschillinel 1 1 977^ . Thus before solid material 
has orbited the planet once, it will fall onto the planet. 
The fastest drifting solids are given by unity in the pa- 
rameter ttpTf = 1, with f2 p the Keplerian frequency and 
Tf the friction time of the solid body. For the physical con- 
ditions in our model DR4 (e.g. density and temperature) 
the fastest drifting objects would be 1 meter sized boul- 
ders, which drift into Jupiter in about a month if released 
at 0.1 Roche radii. Thus the drift problem for the forma- 
tion of moons from circumplanetary material is even more 
crucial than in the case of planet formation. We conclude 
that at least for the time when Jupiter is accreting its 
mass there can be no formation of a satellite system. 

4. Conclusions 

We have performed full 3D radiation hydrodynamical sim- 
ulations of embedded protoplanets in disks, and compare 
the results in detail to the standard isothermal approach. 

Mean torques and migration rates are not strongly af- 
fected by the treatment of the thermodynamics in the case 
of Jupiter mass planets. This might change for lower mass 
planets, which are more embedded in the disk. The fluctu- 
ations of the torques on the other hand are much stronger, 
in particular in higher resolution cases. Similar effects 
have been observed in high resolution nested grid simu- 
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Fig. 15. Model DR4: Temperature in the r-0 plane of the protoplanetary disk at the azimuthal location of the planet 
after 141 orbits. Brightness is logarithmic temperature (lighter = warmer) between 30K and 1500K, contours are 
equi-density lines (in g cm -3 ) and vectors give the logarithmic mass flux. See the online edition for a color version of 
this plot. 
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Fig. 14. Model DR4: Pressure scale height in the mid- 
plane of the circumplanetary cloud after 141 orbits. Radial 
location is given in units of Hill radii. 

lations l|D'Angelo et al.ll2002t l2003bl) and also MHD sim- 



ulations of planet-disk interactions JPapaloizou fe Nelson! 
120031 iNelson fc PapaloizovJl2003l Iwinters et al.ll2003|) . In 

some cases the torques even change their sign for a short 
period. 

We find that planets arc most likely to form a circum- 
planetary pressure-supported envelope rather than an ac- 
cretion disk around them, with strong convective vertical 
flows. The relative pressure scale height in the circumplan- 
etary material is at least 0.5, in which case the approxi- 
mations for a thin Keplerian accretion disk are no longer 
valid. What results is a cloud which rotates at only 50 
percent of the Keplerian value. 

This slow rotation of the circumplanetary gas makes 
it impossible to form the moons of Jupiter from circum- 
planetary material as the solid content will rain down the 
planet in less than an orbit around Jupiter. 

How does this thick and slowly rotating envelope be- 
have with respect to magneto-centr ifugal wind a s dis- 
cussed for a circumplanetary disk by iFendtl fcOOaft ? The 
high thermal pressure could be beneficial for launching 
the magneto-centrifugal wind, but on the other hand we 
observe a strong mass infall onto the poles of the planet 
and the surface of the circumplanetary cloud, thus a wind 
might be suppressed. Detailed investigations on that issue 
are needed. 

The deviation from symmetry with respect to the mid- 
plane is very small in our simulations. Some convection 
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cells show eventual overshooting across the midplane but 
in general, treating only one hemisphere would have led 
to the same results. 

We observe strong vertical flows in the early gap open- 
ing phase. Considering the entire mass accretion phase 
of the young; plane t starting from a few earth masses 
I Pollack et al.l fl9961 . one finds that these vertical foun- 
tains may last for hundred-thousands of years before the 
planet opens its gap, making the effect clearly relevant for 
observations. As a result there would be locally a stronger 
flaring of the disk, and more radiation could interact with 
the small dust grains in the surface layer above the planet. 
Future instruments like SIM should be able to observe 
this asymmetry in the scattered light from disks in which 
planets are forming (Geoffrey Bryden personal communi- 
cation). On the other hand, the effect might be weaker if 
we allow for a gradual buildup of the planet. This has to 
be investigated in future simulations. 

We do not observe the formation of vortices at the 
edges of the gap, which can often be found in th e 
equivalent 2D simulations (|Klevlll999UKoller et afll2003|) . 
Further investigations will have to show whether this is re- 
lated to the 3D nature of our simulations or a consequence 
of the comparable low resolution of the grid. 

Finally, we speculate that the extended warm gas and 
dust which fills most of the Roche lobe should be dc- 
tectable in the far infr ared. Investigations such as by 
IWolf fc D'Aneelol l|2005|) indicate this already for 2D flat 
simulations that include radiative cooling. Similar studies 
using 3D density and temperature distributions will follow 
in the near future. 
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-4 


5.5 


X 


10" 


-4 


2 





X 


io- 


11 


800 


1.0 


X 


10 5 


DR4 


100 


X 


25 


X 


200 


lc 


rad+ 





1 


0.1 




141 


3.5 x 10" 


-4 


4.5 


X 


10" 


-4 


4 


5 


X 


10" 


10 


1500 


1.2 


X 


10 5 



Table 1. Parameters chosen in the different simulations and their results. These parameters are the dimensioning of the grid 
(n r ,n$,n v ), the kind of spacing (lc = logarithmically centered around the planet), the thermodynamics (iso = locally isothermal, 
rad = ideal gas plus radiation transport, rad+ = ideal gas, radiation transport and conservation of total energy for gas accreted 
onto planet), the gravitational smoothing radius (q s = fraction of the Hill radius over that smoothing is applied), the accretion 
radius (g acc = fraction of the Roche lobe (respectively Hill radius) that mass is taken out of), the type of the initial model (G: 
starting with an existing gap), and the total number of orbits calculated (tint). The results we give are the accretion rate onto 
the Planet (M) in units of Jupiter masses per orbit, the mass contained in the Roche lobe (Mrtoche)) the maximum density and 
temperature in the Roche lobe (p m ax,2max) and the migration time-scale that we determined (r m i g ). 



